/*
 * Membrane.cpp
 *
 *  Created on: Nov 9, 2011
 *      Author: yi
 */

#include "mbuild/Membrane.h"
const double Membrane::Inf=1e298;

Membrane::Membrane() {
	myfunc=NULL;
	funcdensitydist=NULL;
	holefunc=NULL;
}

Membrane::~Membrane() {
	delete holefunc;
}

void Membrane::set_cutoff(double cutoff)
{
	this->cutoff=cutoff;
}



void Membrane::set_distribution(Dfunc *f,InhomogenuousFunc* d)
{
	/*
	 * set Distribution function & Density distribution.
	 */
	myfunc=f;
	funcdensitydist=d;
}



Dfunc *Membrane::get_distribution_fun(double x, double y, double z)
{
	return myfunc;
}



double Membrane::density(double x0, double y0, double z0, double dx, double dy, double dz)
{
	/* Calculate the integration of density defined by myfunc* in a grid box.
	 * This grid box is centered at x0,y0,z0
	 * and has width of dx,dy,dz on x,y,z axis
	 *
	 */
	double pos[MAX_POINTS][3];
	int npoint=0;
	double a,b,c,d;
	double x,y,z,dist;
	int ncyc=3;
	if (myfunc == NULL) throw NoAttachedDistributionFunc();
	closet_points(x0,y0,z0,pos,npoint); //find the closet point to (x0,y0,z0) on membrane surface
//	PRINT("%5.2f %5.2f %5.2f\n",x0,y0,z0);
//	PMAP(pos,1,3,"%5.2f");
	npoint=1; //For current implementation, just get the first point.
	double sum=0;
	for(register int i=0;i<npoint;i++) {
		x=pos[i][0];
		y=pos[i][1]; 
		z=pos[i][2];
		dist=sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0));
		if (in_rim(x,y,z) && ( dist > 0 )) {
		 //for discontinuous curved surface, the normal of tangent plane in rim is (x,y,z)-(x0,y0,z0)
          a=x0-x;b=y0-y;c=z0-z;
		} else {
		  a=df_dx(x,y,z);b=df_dy(x,y,z);c=df_dz(x,y,z);
		}
		d=-(a*x+b*y+c*z); // calculate the tangent plane at closet point (x,y,z)
		myfunc->set_plane(a,b,c,d); // approximate the curved surface into its tangent plane
		if (funcdensitydist != NULL) {
			funcdensitydist->modify(myfunc,x,y,z);
		}
		if (dist<= cutoff) {
			// only when distance is less than cutoff, the density is the average over the whole volume.
		  sum+=integrate(*myfunc,x0-dx/2,x0+dx/2,y0-dy/2,y0+dy/2,z0-dz/2,z0+dz/2,ncyc,ncyc,ncyc)/(dx*dy*dz);
		} else {
			// otherwise the density is just the same as the center of the grid.
		  sum+=(myfunc->f(x0,y0,z0));
		}
	}
	// If the calculating property is charge density or potential density(which are addictive), the density is the sum of contribution from several neighboring surface.
	if (!(myfunc->IsAdditive) && npoint>0) sum/=npoint; //If the density is dielectric constant which is not addictive, the density is the average of contribution from several neighboring surface.
//	PVAR(sum);
	return sum;
}



void Membrane::set_integration_cyc(int ncyc)
{
	this->ncyc=ncyc;
}



void Membrane::set_region(double xl, double xh, double yl, double yh, double zl, double zh,bool dofastmatch,double dx,double dy,double dz)
{
	/*
     * Set the region inside which we build the membrane.
     * If dofastmatch is true, fast matching method is used to calculate the distance of points to membrane surface.
	 */
	this->xl=xl;
	this->xh=xh;
	this->xcen=(xl+xh)/2;
	this->yl=yl;
	this->yh=yh;
	this->ycen=(yl+yh)/2;
	this->zl=zl;
	this->zh=zh;
	this->zcen=(zl+zh)/2;
	if (dofastmatch) { fast_match(dx,dy,dz);}
}



bool Membrane::in_region(double & x, double & y, double & z)
{
	  // check if point (x,y,z) is in the building region.
	  if (x>=xl&& x<=xh && y>=yl && y<=yh && z>=zl && z<=zh) {
	    x-=xcen; //shift x,y,z to the local coordinates of the building region;
	    y-=ycen;
	    z-=zcen;
	    return (true);
	  } else return (false);
}



GridSpace *Membrane::fill(GridSpace *tofill, GridSpace *reference)
{
	/* Fill the GridSpace *tofill with defined Membrane
	 * a reference map is used to help identify the grid points that has already been occupied by solute atoms.
	 */
	unsigned long nx,ny,nz,nx2,ny2,nz2,n;
	n=tofill->get_dimension(nx,ny,nz);
	reference->get_dimension(nx2,ny2,nz2);
	if (nx!=nx2 || ny!=ny2 || nz!=nz2) throw GridSpaceNotEqlRef();
	double dx,dy,dz,x0,y0,z0;
	tofill->get_grid_size(dx,dy,dz);
	tofill->get_origin(x0,y0,z0);
	double* ref=reference->values();
	//#pragma omp parallel for
	for(int ix=0;ix<nx;ix++)
	for(int iy=0;iy<ny;iy++)
	for(int iz=0;iz<nz;iz++)
	{
		unsigned long  i=ix*ny*nz+iy*nz+iz;
		double x=ix*dx+x0;
		double y=iy*dy+y0;
		double z=iz*dz+z0;
		// only fill grids inside build region
		if (in_region(x,y,z)){
			if(myfunc->IfEqualSolute(ref[i])) continue; //if grid already have value equal to skip value of myfunc, then do nothing. This is for skipping grids that already contained solute data.
			else if (! myfunc->IsAdditive)  (tofill->data)[i]=((this->density)(x,y,z,dx,dy,dz)); // for dielectric constant which is not additive just set the new value.
			else (tofill->data)[i]+=((this->density)(x,y,z,dx,dy,dz)); //for charge density... which is additive, add the new density
		};
//	  PRINT("%5.2f %5.2f %5.2f %8e\n",x,y,z,tofill->data[i]);
	}
	return (tofill);
}



void Membrane::fast_match(double dx, double dy, double dz)
{
	/*
	 *  Build a distance matrix inside building region using fast matching method
	 *  The building region is discretized  into grids with width dx*dy*dz.
	 */
	double x,y,z,x1,y1,z1,d,min;
	unsigned long nx,ny,nz;
	register unsigned long i,newi,j,t,imin,sx,sy,sz,ixx,iyy,izz;
	const int R=4; //set search radius on the surface.
	// set up the grid space to store distance and source point(X0,Y0,Z0)
	// it will be slightly larger than xh,yh,zh to ensure interpolation can be used to get value.
	nx=1+ceil((xh-xl)/dx);ny=1+ceil((yh-yl)/dy);nz=1+ceil((zh-zl)/dz);
	distance.set_dimension(nx,ny,nz);
	distance.set_origin(xl,yl,zl);
	distance.set_grid_size(dx,dy,dz);
	X0.set_dimension(nx,ny,nz);Y0.set_dimension(nx,ny,nz);Z0.set_dimension(nx,ny,nz);
	X0.set_origin(xl,yl,zl);Y0.set_origin(xl,yl,zl);Z0.set_origin(xl,yl,zl);
	X0.set_grid_size(dx,dy,dz);Y0.set_grid_size(dx,dy,dz);Z0.set_grid_size(dx,dy,dz);
	distance.initialize(Inf);X0.initialize(Inf);Y0.initialize(Inf);Z0.initialize(Inf);

	queue<unsigned long> positiveband;
	queue<unsigned long> negativeband;
	queue<unsigned long> zeroband;
	unsigned long* source;
	try {
		source=new unsigned long[nx*ny*nz];
	} catch(bad_alloc) {
		throw NotEnoughMemory();
	}
	if (source == NULL ) throw NotEnoughMemory();
	//Guess the zero band
	for(register int ix=0;ix<nx;ix++) {
		for(register int iy=0;iy<ny;iy++) {
			for(register int iz=0;iz<nz;iz++) {
				x=ix*dx+xl;y=iy*dy+zl;z=iz*dz+zl;
				i=ix*ny*nz+iy*nz+iz;
				// If the point(x,y,z) is on the surface(or close enough to the surface).
				if (fabs(f(x,y,z))<fabs((df_dx(x,y,z)*dx/2)+fabs(df_dy(x,y,z)*dy/2)+fabs(df_dz(x,y,z)*dz/2))) {
					X0.data[i]=x;Y0.data[i]=y;Z0.data[i]=z;
					source[i]=i;
					distance.data[i]=0.0;
					zeroband.push(i);positiveband.push(i);negativeband.push(i);
					//printf("%lf %lf %lf\n",x,y,z);
				}
				if (f(x,y,z)<0) { distance.data[i]=-Inf;}
			}
		}
	}
	//Matching towards positive distances
	while (positiveband.size() > 0) {
		//remove grid with known distance from narrow band.
		i=positiveband.front();
		positiveband.pop();
		//closet neighbors of i
		for(register int ix=-1;ix<=1;ix++)
			for(register int iy=-1;iy<=1;iy++)
				for(register int iz=-1;iz<=1;iz++)
				{
					newi=i+ny*nz*ix+nz*iy+iz;
					if(i == newi) continue;
					ixx=ix+i/(ny*nz);iyy=iy+(i/nz)%ny;izz=iz+i%nz;
					if( ixx <0 || ixx >= nx || iyy<0 || iyy >=ny || izz <0 || izz >=nz) continue;
					x=ixx*dx+xl;
					y=iyy*dy+yl;
					z=izz*dz+zl;
					if ( distance.data[newi] == Inf ) {
						//push undefined datapoint into narrow band.
						positiveband.push(newi);
					} else {
						continue;
					}
					//printf("x %lf y %lf z %lf v %lf distance %e\n",x,y,z,f(x,y,z),distance.data[newi]);
					if(distance.data[newi] > distance.data[i]) {
						//Recomputing the distance based on i's neighbor
						min=Inf;
						imin=0;
						sx=source[i]/(ny*nz);
						sy=(source[i]/nz)%ny;
						sz=source[i]%nz;
						for(register int jx=-R;jx<=R;jx++)
							for(register int jy=-R;jy<=R;jy++)
								for(register int jz=-R;jz<=R;jz++){
									if (sx+jx>=0 && sx+jx<nx &&	sy+jy>=0 && sy+jy<ny && sz+jz>=0 && sz+jz<nz) {
										j=source[i]+jx*ny*nz+jy*nz+jz;
										if(distance.data[j]!=0.0) continue;
										x1=(j/(nz*ny))*dx+xl;
										y1=((j/nz)%ny)*dy+yl;
										z1=(j%nz)*dz+zl;
										d=(x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1);
										if (d< min) {min=d;imin=j;}
									}
								}
						//printf("min %lf\n",min);
						//define distance of neighbor j
						distance.data[newi]=sqrt(min);
						source[newi]=imin;
						X0.data[newi]=(imin/(nz*ny))*dx+xl;
						Y0.data[newi]=((imin/nz)%ny)*dy+yl;
						Z0.data[newi]=(imin%nz)*dz+zl;
					}
				}
	}
	//Matching towards negative distances
	while (negativeband.size() > 0) {
		//remove grid with known distance from narrow band.
		i=negativeband.front();
//		printf("Pop %i\n",(int)i);
		negativeband.pop();
		//closet neighbors of i
		for(register int ix=-1;ix<=1;ix++)
			for(register int iy=-1;iy<=1;iy++)
				for(register int iz=-1;iz<=1;iz++)
				{
					newi=i+ny*nz*ix+nz*iy+iz;
					if(i == newi) continue;
					ixx=ix+i/(ny*nz);iyy=iy+(i/nz)%ny;izz=iz+i%nz;
					if( ixx <0 || ixx >= nx || iyy<0 || iyy >=ny || izz <0 || izz >=nz) continue;
					x=ixx*dx+xl;
					y=iyy*dy+yl;
					z=izz*dz+zl;
					if ( distance.data[newi] == -Inf ) {
						//push undefined datapoint into narrow band.
						negativeband.push(newi);
					} else {
						continue;
					}
					//printf("x %lf y %lf z %lf v %lf distance %e\n",x,y,z,f(x,y,z),distance.data[newi]);
					if(distance.data[newi] < distance.data[i]) {
						//Recomputing the distance based on i's neighbor
						min=Inf;
						imin=0;
						sx=source[i]/(ny*nz);
						sy=(source[i]/nz)%ny;
						sz=source[i]%nz;
						for(register int jx=-R;jx<=R;jx++)
							for(register int jy=-R;jy<=R;jy++)
								for(register int jz=-R;jz<=R;jz++){
									if (sx+jx>=0 && sx+jx<nx &&	sy+jy>=0 && sy+jy<ny && sz+jz>=0 && sz+jz<nz) {
										j=source[i]+jx*ny*nz+jy*nz+jz;
										if(distance.data[j]!=0.0) continue;
										x1=(j/(nz*ny))*dx+xl;
										y1=((j/nz)%ny)*dy+yl;
										z1=(j%nz)*dz+zl;
										d=(x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1);
										if (d< min) {min=d;imin=j;}
									}
								}
						distance.data[newi]=-sqrt(min);
						source[newi]=imin;
						X0.data[newi]=(imin/(nz*ny))*dx+xl;
						Y0.data[newi]=((imin/nz)%ny)*dy+yl;
						Z0.data[newi]=(imin%nz)*dz+zl;
					}
				}
	}
	//DxFile::write_dx("distance.dx",&distance);
	delete source;
}

inline double Membrane::fast_sign(double f) {
	/*
	 * Just to get the sign of a float number.
	 */
    if (f > 0) return 1;
    return (f == 0) ? 0 : -1;
}

void Membrane::closet_points(double x, double y, double z, double points[][3], int & npoint)
{
	/*
	 *  Default subroutine: get the closet_points from pre-calulated matrix obtained from fast_match method.
	 */
	points[0][0]=X0.getValue(x,y,z);
	points[0][1]=Y0.getValue(x,y,z);
	points[0][2]=Z0.getValue(x,y,z);
	//printf("%f %f %f\n",points[0][0],points[0][1],distance.getValue(x,y,z));
	npoint=1;
}
void Membrane::has_hole(HoleFunc* f) {
	//TODO Notice: shallow copy of object. If the object is declared inside different scope. Will give segment fault if use it in another scope.
	this->holefunc=f->clone();
}

GridSpace *Membrane::create_hole(GridSpace *tocreate, GridSpace *reference)
{
	/*
	 * Create a hole inside a region(where holefuc gives negative number) of GridSpace.
	 * Inside the hole, the grid is filled with density values at Infinity distance from surface.
	 * In the case of membrane, it is the density values(charge,dielectric constant...etc.) of water.
	 * The grid points that already have been occupied by solute atoms will be skipped.
	 * A reference grid space is used to tell which grid point is occupied by solute atoms.
	 */
	unsigned long nx,ny,nz,rnx,rny,rnz,i;
	double dx,dy,dz,x,y,z,x0,y0,z0,xh,yh,zh;
	double rx0,ry0,rz0,rdx,rdy,rdz,rxh,ryh,rzh;
	// Improve? store xh,yh,xl,yl...etc in GridSpace? or provide a function in GridSpace?
  	unsigned long n=tocreate->get_dimension(nx,ny,nz);
  	tocreate->get_grid_size(dx,dy,dz);
  	tocreate->get_origin(x0,y0,z0);
    reference->get_dimension(rnx,rny,rnz);
  	reference->get_grid_size(rdx,rdy,rdz);
  	reference->get_origin(rx0,ry0,rz0);
  	xh=nx*dx+x0;yh=ny*dy+y0;zh=nz*dz+z0;
  	rxh=rnx*rdx+rx0;ryh=rny*rdy+ry0;rzh=rnz*rdz+rz0;
  	// skip if no function is given for the hole
  	if (holefunc == NULL) return tocreate;
  	// reference GridSpace must equal to the computing GridSpace. If not, raise an error.
	if (nx!=rnx || ny!=rny || nz!=rnz) throw GridSpaceNotEqlRef();
	for( unsigned int ix=0;ix<nx;ix++)
	for( unsigned int iy=0;iy<ny;iy++)
	for( unsigned int iz=0;iz<nz;iz++)
	{
		// Option: store pre-calculated x,y,z in grid space will be faster?
		//	  x=tofill->x[i];
		//	  y=tofill->y[i];
		//	  z=tofill->z[i];
		i=ix*ny*nz+iy*nz+iz;
		x=ix*dx+x0;
		y=iy*dy+y0;
		z=iz*dz+z0;
//		printf("%f\n",holefunc->f(x,y,z));
		if (holefunc->f(x,y,z) < 0) {
			if(!(myfunc->IfEqualSolute(reference->data[i] ))) {tocreate->data[i]=myfunc->f(Inf);}
		}
	 }
	 return tocreate;
}



